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Transverse Mercator with an accuracy of a few nanometers 
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Implementations of two algorithms for the transverse Mercator projection are described; these achieve accuracies 
close to machine precision. One is based on the exact equations of Thompson and Lee and the other uses an 
extension of Kriiger's series for the mapping to higher order. The exact method provides an accuracy of 9 nm 
over the entire ellipsoid, while the errors in the series method are less than 5 nm within 3900 km of the central 
meridian. In each case, the meridian convergence and scale are also computed with similar accuracy. The speed 
of the series method is competitive with other less accurate algorithms and the exact method is about 5 times 
slower. 
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1. INTRODUCTION 

The transverse Mercator or Gauss-Kriiger projection is a 
conformal mapping of the earth ellipsoid where a central 
meridian is mapped into a straight line at constant scale. Be- 
cause it cannot be expressed in terms of elementary functions, 
the mapping is usually computed by means of a truncated se- 
ries (Krugeii lTQll I^ oma si Il952f) . The resulting mapping 
approximates the true mapping only within a region centered 
on the central meridian. 

Transverse Mercator is one of the commonest projections 
used for large-scale maps (it is used for the grid systems 
of several countries and is the basis of the univ ersal trans- 
verse Mercator (UTM) system dHager ef"aZ I EM, Chap. 2)). 
For the WGS84 ellipsoid, the variation of the scale is 1.25% 
within 1000 km of the central meridian; it is therefore desir- 
able to find algorithms for the mapping which are accurate 
to machine precision over at least this area. In this paper, 
I describe the implementation of two such a lgorith ms, one 
based on the exact equations given by L ed (Il976h and the 
other extending the series given by IKriigeii (Il912h to higher 
order. Both implementations compute the forward and re- 
verse mappings and also return the meridian convergence and 
scale. Th ese implementations are included in GeographicLib 
(iKarnevl 12010) . 

Scores of other authors have presented methods for comput- 
ing this mapping over the past century. In particular, iDoziej 
([l980) provided an implemen tation of Lee's exact method, 
and lEngsager and Podeil (l2007l) give Kriiger's series to 7th or- 
der. The distinguishing aspects of this work are the reduction 
of the overall numerical errors (truncation and round-off) to 
close to the precision limit of the computer and the concrete 
bounds 1 place on these errors. 
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Because floatin g-point numbers have a finite spacing 
(lOlver et al. 1 12010'. p.l(i)), the limiting accuracy of any im- 
plementation is about AI/2P where M — 10 000 km is the 
length of the quarter meridian of the earth and p is the num- 
ber of bits in the fraction of the floating-point number system. 
This gives an error limit of 0.5 m forp = 24 (single precision 
or float), 1 nm for p = 53 (double precision or double), and 
0.5pmforp = 64 (extended precision or long double). (Here, 
I use SI prefixes: 1 nm = 10~^ m, 1 pm = 10~^^ m.) Typi- 
cally p = 24 is too inaccurate to be useful and 1 don't consider 
this further in this paper. My standard working precision is 
double and the resulting accuracy, if it can be achieved, would 
satisfy most needs. However, 1 also use extended precision as 
one of the tools to verify the accuracy of the double precision 
implementations . 

Formulas for mappings can contain expressions which are 
numerically ill-conditioned causing precision to be lost. This 
loss of precision is of little consequence if the truncation er- 
rors are of the same order However, in attempting to min- 
imize the numerical errors, 1 needed an accurate means of 
quantifying the truncation and round-off errors. To this end, 
1 constructed a large test set of projected points which were 
computed with an accuracy of 80 decimal digits. This allowed 
me to eliminate many sources of round-off error. The resulting 
accuracies are about 4-8 times the limiting value (equivalent 
to a loss of only 2-3 bits of precision) and this applies to both 
double and extended precisions. 

In Sect. 121 I review the series method given by iKriigeil 
(I1912I) modifying it to minimize the round-off errors. I turn 
next. Sect. [3] to the formulation of the exact transverse Mer- 
cator projection bv iLeel ([1976) which I use to construct the 
high-precision test set; I also describe its implementation us- 
ing double precision and I quantify the round-off errors. I 
extend Kriiger's series to 8th order (see Sect. |4| and give the 
truncation error for the series as a function of truncation level 
and distance from the central meridian. Finally, in Sect. [Sj I 
discuss some of the properties of the exact mapping far from 
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the central meridian. 



2. KRUGER'S SERIES 



I summarize here the method developed by iKrugerl (Il912l 
§§5-8), simplifying it and adapting it for optimal implemen- 
tati on on a computer. The method is also briefly described 
by iBugavevskiv and Snvde3 (Il995' , ii5.1.6). The method en- 
tails mapping the ellipsoid to the conformal sphere and for 
this reason I begin by describing the spherical transverse Mer- 
cator projection. 

Consider a sphere and a point on that sphere of latitude 
(f)' and longitude relative to the central meridian of A. (I use 
primes on variables, e.g., </)', where necessary, to distinguish 
them from their ellipsoidal counterparts.) The isometric lati- 
tude is given by 



(1) 



sech tdt = tan sinh x = sin tanh : 



where 

gdx = 

Jo 

is the Gudermannian function given by lOlver et al.\ ilQld 
f:4.23(viii)) (henceforth referred to as DLMH , the Digital Li- 
brary of Mathematical Functions) and 



gd 



sec tdt — sinh ^ tan x — tanh ^ sin x 



is its inverse. The standard (equatorial) Mercator projection 
maps the sphere onto the plane (A, When working with 
conformal mappings it is often useful to represent coordinates 
with complex numbers where the real part represents the nor- 
thing and the imaginary part the easting (the phase, or argu- 
ment, of the complex number gives a bearing measured clock- 
wise). In this representation the Mercator projection (a con- 
formal mapping) is given by 

Any analytic function of x ^Iso represents a conformal map- 
ping (except where its derivative vanishes); its derivative gives 
the change in the meridian convergence and scale for the map- 
ping. In particular (.Lee, ,1976 , Eq. (12.3)), 



(2) 



gives the transverse Mercator projection of the sphere. This 
is easy to confirm by evaluating the mapping for A = 0; this 
gives (' — (/)', i.e., the central meridian is mapped to a straight 
line at constant scale (the defining property of the mapping). 

I consider now an ellipsoid of revolution with equatorial 
radius a, polar semi-axis b, flattening / = {a — b)/a, eccen- 
tricity e = ^ /(2 — /), and third flattening n = {a — b)/ 
{a + b) = //(2 — /). For a point wit h latitude (p and longitude 



A, the isometric latitude is given by (lLambertl[l772[ §117) 



ip = logtan( - + - 



1 , / 1 + e sin ( 
-elog — 

2 V 1 — e sm t 



Using the identities 'dlmf", Eqs. |(4.23 .42)1 and [(43T24)1 this 
relation may also be written as 



tp ~ gd ^ (f) — e tanh ^ (e sin ( 



(3) 



As in the case of the sphere, x = i/) + iA defines the Mercator 
projection. Equating the isometric latitude for the sphere with 
that for the ellipsoid, ijj' — ij}, defines a relation 



y = gd(gd "^0— etanh "^(esin^)), 



(4) 



which maps a point on the ellipsoid with latitude </> confor- 
mally to a point on the sphere with latitude 0'. In this con- 
text, (f)' is called the "conformal latitude" and the sphere is 
referred to as the "conformal sphere." The transformation to 
C', Eq. (|2]l, where ip is given by Eq. ([3]), defines a confor- 
mal mapping of the ellipsoid to a plane in which the central 
meridian is mapped to a straight line with a scale which is 
nearly constant (the variation is 0(f)). I call this the "spher- 
ical transverse Mercator projection" as it is the simplest gen- 
eralization of the spherical projection to the elhpsoid. The 
graticule for this mapping is shown in Fig.[TTa). lKrugeii ( ll912i 
§8) now "rectifies" this mapping by applying a near-identity 
transformation to to make the scale constant which yields 
the Gauss-Ki-iiger mapping ( — irj with 



c 



C' + ^ajSin2jC', 



(5) 



where aj is real (this form for the transformation is derived 
in Sect. |4|l. Similarly the transformation from C,' to ^ can be 
written as 



c' = c- 



oo 



(6) 



where I3j is real. Rriiger's expressions for aj and p)j are given 
below, Eqs. (fT2b and ( fTTI l. and we outline their derivation in 
Sect.m 

First, I address the computation of given (j> and A with 
special emphasis on maintaining numerical accuracy. Follow- 
ing Kriiger, I write C' = + W and give separate equations 
for ^' and rj'. However, in order to maintain accuracy near 
(j> — zL^TT, I use T — tani/) and r' = tant/)' and eliminate 
(f)' and ijj from the relations. An expression for r' is found by 
taking the tangent of Eq. (|4|l and using the addition rule for 
the hyperbolic sine to give 



r' = Tvr 



where 



tan ( 



(7) 

(8) 
(9) 



a — sinh (e tanh ^{er/ \/\ + t^) ) . 

Eliminating p' from the expressions for ^ and vl dRriigeri 
il91Z Eq. (8.36)) yields 

£' = tan" ^(t7 cos A), 

(10) 



77' = sinh ^ (sin A/ \/r'2 + cos^ A) . 
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Splitting Eq. (|5]l into real and imaginary parts gives 

oo 

e = e' + ^aj sin2je' cosh2jV, 

oo 

T] — rj' + Oj cos 2j^' sinh 2jr/', 

where (lKrugeilll912[ Eq. (8.41)) 

1 2 2 5 3 41 4 

ai — —n n H n H n + • • • , 

2 3 16 180 

13 9 3 Q 55T ^ 



as 



61 o 103 4 
-n —n 



240 
49561 

161280 ' 



140 



(11) 



(12) 



Finally, ^ and rj are scaled to give the transverse Mercator 
easting x and northing y. 



Equation (|7]l may be inverted by Newton's method, 

for i = 0, 



ri_i + (3T,;_i, Otherwise, 



t[ = ryi + o-f- cryi + rf, 



VT+Tf (1 -e2)VrT 



(20) 
(21) 



This usually converges to round-off after two iterations, i.e.. 
T = T2, which gives 



tan T. 



(22) 



The meridian convergence and scale can be found during 
the forward mapping by differentiating Eq. ^ and writing 



(13) 



where fcg is the scale on the ce ntral meridian, 2'kA is the cir- 
cumference of a meridian, and (lKrugeiill912l Ea. (5.5)) 



1+71 



+64" 



(14) 



p' = 1 + ^ 2jaj cos 2j^' cosh 2j?7', 

oo 

g' = 2jaj sin 2j^' sinh 2jr7'. 



(23) 



Typically fco is chosen to be slightly less than 1 to minimize 
the deviation of the scale from unity in some region around 
the central meridian. 

Converting from transverse Mercator to geographic coordi- 
nates entails reversing these steps. Equations (flJl i give 

T^^x/(k„A), ^ = y/ikoA). (15) 

lKrugeil(ll912[ §7) writes (' in terms of ( by inverting Eq. (fTTl i 
to give 



Then the meridian convergence (the bearing of grid north, the 
y axis, measure d clockwise fro m true north) is given by 7 = 
7' + 7", where (lKrugeiill912[ Eqs. (8.44^5)) 



7' = tan-i((r7Vl + r'^) tanA), 
7" = tan-i(g7p'). 



(24) 



The scale is given by fc = k^k'k" , where jKriigeil Il912l 
Eq. (8.47)) 



^' = ^ - ^ /3j sin cosh 2j-q, 

00 

■q' — fj — j3j cos 2j^ sinh 2jrj, 



where (lKrugeilll912[ Eq. (7.26*)) 



n 1 S 1 4 

Pi = — n n H n n 

^ 2 3 96 360 



1 



1 



437 



= + 15"' - 1440' 



480 840 



Pi 



180 
4397 



161280 

Inverting Eq. gives dKrugeii [19121 Eq. (7.25)) 



t' sin V sinh^ 77' + cos^ 
A = tan^^(sinhry'/ cos^'). 



(16) 



(17) 



(18) 



k' = VT — e" sm 



2 „;„2 ^^l + r2/v/T'2 + C0s2A, 

A , (25) 

k" = -Vp'^ +q'2. 
a 

Here 7' and fc' give the convergence and scale for the spheri- 
cal transverse Mercator projection, while 7" and fc" give the 
corrections due to Eqs. (|5]l and ( fT3b . 

To determine the convergence and scale during the reverse 
mapping, differentiate Eq. ^ and write 



p + iq 



dC 1 



or 



p = 1 — 2j/3j cos 2j^ cosh 2j?7, 

00 



(26) 



4 



The c onvergence is given by 7 = 7' + 7", where dKriigeii 
[19121 Eqs. (7.31-31*)) 



7' = tan ^(tan^' tanh?7'), 
7" = tan-i(g/p). 



(27) 



The scale is given by fc = kok'k", where dKriigeii Il912l 
Eq. (7.33)) 



k' = Vl-e'^ sin^ 0^1 + sinh^ t]' + cos^ C', 

A 1 (28) 



fc" 



In summary, Kriiger's methods for the forward and re- 
verse mappings are given by the numbered Eqs. d7ll-(fl4li and 
Eqs. (fl4]l-d22]l. respectively. The scale and meridian conver- 
gence are similarly given by the Eqs. (|23]l-(|25]l during the for- 
w ard mapp ing and Eqs. (l26ll-d28Tl during the reverse mapping. 

iKrii ger' truncates the series at order n*, as shown here. This 
results in very small errors, considering that Kriiger published 
his paper in 1912. The maximum of the errors for the forward 
and reverse mappings (both expressed as true distances) is 
0.31 /im within 1000 km of the central meridian and is 1 mm 
within 6000 km of the central meridian. The truncated map- 
ping is exactly conformal; however Eqs. ^ and ^ are not 
inverses of one another if the sums are truncated. It is, of 
course, possible to construct an exact inverse of the truncation 
of Eq. (|5]l, e.g., by solving it using Newton's method. How- 
ever, in practice, it is better merely to retain enough terms in 
the sum so that the truncation error is less than the round-off 



error. 



In numerically implementing this method, the terms A, aj, 
and 13 j, Eqs. ( fT4] i. (fT2] i. and dlTl i. need only be computed once 
for a given ellipsoid and, for acc uracy and speed, should be 
evaluated in Horner form (IdlmfI §1.11(1)) ); for example, ai 
is evaluated to order n"^ as 



ai 



+ (-2- + ALn) 

^ VI6 ^ 180 



Furthermore the trigonometric series, Eqs . ( fTTT i. ( fTSI l. ( l23T l, 
and (|26l l. can be evaluated using IClenshawl(ri955) summation 
(IdlmfI i3.11(ii)| l which minimizes the number of evaluations 
of trigonometric and hyperbolic functions. Thus Eqs. ( fTTT l and 
may be summed to order J with 

Cj+l = C,7+2 = 0, 

Cj = 2cj+i cos2(^' + irj') - Cj+2 + aj, 
£. + i-q ^ i' + irj' + ci siii2(^' + iry'), 



and 



dj+i — dj+2 — 0, 

dj = 2dj+i cos2(^' + irj') — 
' + iq' = l~d2+di cos 2(^' + irj'), 



separated into real and imaginary parts and with the recursion 
relations for Cj and dj evaluated for J > j > 0. The summa- 
tions of Eqs. (fTST i and d26T l are handled in a similar fashion. 



My introduction to Kriiger's expansion was a rep ort by 
the Finnish Geodetic Institute dKuittinen et al. I l2006b . The 
method described here follows this report with a few changes 
to improve the numerical accuracy: (a) I use more stable for- 
mulas for converting from geographic to the spherical trans- 
verse Mercator coordinates; (b) I solve for the geographic lati- 
tude by Newton's method instead of by iteration; and (c) I use 
Kriiger's method for determining the convergence and scale 
instead of less accurate expansions in the longitude. 

In contrast to the series given here, the formulas given by 
Kriiger in a later section of his paper, §14, involve an expan- 
sion in the longitude difference instead of the flattening. This 
expansion forms the basis of the approximate transverse Mer- 
cator formulas pres ented by ThomasI (1 19521 pp. 2-6) and in 
th e report on UTM (jHager et ai, 1989, Chap. 2) and are used 
in lGeotransI d2010l) . For computing UTM coordinates, the er- 
rors are less than 1 mm. Unfortunately, the truncated series 
does not define an exact conformal mapping. In addition, in 
some applications, use of these series may lead to unaccept- 
ably large errors. For example, consider mapping Greenland 
with transverse Mercator with a central meridian of 42° W. 
The landmass of Greenland lies within 750 km of this cen- 
tral meridian and the maximum variation in the scale of trans- 
verse Mercator is only 0.7% — in other words, the transverse 
Mercator projection is ideal for this application. The error in 
computing transverse Mercator with Kriiger's 4th order series 
is (as we have seen) less than 1 fim. However the maximum 
error using Thomas' series (as implemented in Geotrans, ver- 
sion 3.0) is over 1 km. 



3. EXACT MAPPING 

The definition of the transverse Mercator projection given 
at the beginning of Sect.[T]serves to specify the mapping com- 
pletely. (There are two minor qualifications to this statement: 
the central scale, the origin, and the orientation of the cen- 
tral meridian need to be specified; in addition, the mapping 
becomes multi-valued very far from the central meridian as 
detailed in Sect. |5]) Provided that the series in Eqs. (|5]l, dSll, 
and (fl4] i are convergent, the Kriiger series method converges 
to the exact Gauss-Kriiger projection and the truncated series 
are a useful basis for numerical approximations to the map- 
ping. 

There is no problem with the convergence of expression for 
A, Eq. dT4l i. This can be written in closed form as 



A 



2a 



E(e) 



2a 



E{4n/{l + n)^), 



where E{k) is t he comp lete elliptic integral of the second kind 



with modulus k 
in a series using 



DLMF . Eq. (19.2.8)), which may be expanded 
DLMFi Eq.,(19.5.2) to give 



A 



1 



1- 



256" 



25 .n^- 



16384' 



(29) 



This series converges for |n| < 1 and, for small n, the relative 
error in truncating the series is given by the first dropped term. 
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The convergence of Eqs. Q and ^ is more compHcated 
because the sine terms in the summands become large for 
large ry or rf . Indeed, the transverse Mercator projection 
has a singularity in its second derivative at 4> — 0° and 
A — ±Ao where Aq = (1 — e)90° beyond which the series 
will diverge; these points are branch points of the mapping 
(|Whittaker and WatsonL Il927[ §5.7) and the properties of the 
mapping in their vicinity are explored in Sect. |5] In order to 
determine the error in the truncated series , I implem ent the 
formulas for the exact mapping as given by lLeel(ll976 l §§54- 
55) who credits E. H. Thompson (1945) for their development. 
A referee has pointed out t o me that a simi lar formulation was 
independently provided bv lLudwigl (1 1943b . Here 1 give only a 
brief description of Lee's method, referring the reader to the 
document ation and source code for GeographicLib for more 
details (Karnevt 12010) . 

The exact mapping is expressed in terms of an intermediate 
mapping, the Thompson projection, denoted hy w ^ u + iv 
with 0^ 1976, Eqs. (54.5) and (55.5)) 



X 

C 



tanh snw — etanh (esniu), 



(30) 
(31) 



w here sn u is one of the Jacobi elliptic functions with modulus 
e ( IdlmfI §22.2), K{k) is the com plete elliptic integral of the 
first kind with modulus k (IdlmfL Eq. (19.2.8)), and £{x, k) is 
Jacobi's epsilon function (IdlmfI Eq. (22.16.20)). 

When implementing these equations, 1 follow Lee and 
break the formulas in terms of their real and imaginary parts. 
This enables the algorithm to be implemented with real arith- 
metic which allows the expressions to be optimized to mini- 
mize the round-off er ror. The nec essary formulas for Eqs. (|30] | 
and (EB are given by 11163 ( 1197 61 Eqs. (54.17) and (55.4)). 

The computation of the forward (resp. reverse) mapping re- 
quires the inversion of Eq. (|30] | (resp. Eq. dSTTi). 1 perform 
these inversions using Newton's meth od in the co mplex plane. 
The needed derivative o f y i s given bv lLeel(ll976i Eq. (54.21)) 
and d(/dw is given by Lee (1976', Eq. ( 55.9)) which may be 
split into real and imaginary parts with Idlmh Eq. (22.8.3) 
^22.6(iv)| The starting guesses for Newton's method are ob- 
tained by finding approximate solutions using one of three 
methods: (a) by using the limit e ^ 0, (b) by expanding 
about the branch point on the equator (the bottom right corner 
of Fig. |3jc)), or (c) by expanding about the singularity at the 
south pole (the top right corner of Fig. Oc)). (The latter two 
methods require a knowledge of the properties of the mapping 
far from the central meridian; see Sect. |5] for more informa- 
tion.) The most time-consuming task in this implementation 
was optimizing the choice of starting point to ensure that the 
method converges in a few iterations. 1 refer the reader to the 
code for de tails. I also compute the meridian convergence and 
scale usinglLeS (11976*. Eqs. (55.12-13)). 

In order to reduce the round-off errors, I needed to identify 
terms in the formulas with the potential for a loss of precision 
and apply identities for the Jacobi elliptic functions (IdlmfI 
Eq. (22.2.10), §22.6(i)) to recast the formulas into equivalent 
ones with better numerical prop erties. 1 use the procedure 
sncndn given bv lBulirschI (1 1 965h for the elliptic functions and 



algorithms Rp, Rd, and Rq of lCarlsonl(ll995 b for the elliptic 
integrals; these algorithms can yield results to arbitrary preci- 
sion. 

1 provide two implementations of the exact mapping: (a) a 
C++ version using stan dard floa t ing-po int arithmetic and 
(b) an implementation in iMaximal (l2009l) . The latter imple- 
mentation makes use of Maxima's "bigfloat" package which 
permits the calculation to be carried out to an arbitrary pre- 
cision. This was used to construct a large test set for the 
mapping which served to benchmark the C++ implementa- 
tion. This set includes randomly distributed points together 
with additional points chosen close to the pole and other pos- 
sibly problematic points and lines. The mapping is computed 
to an accuracy of 80 decimal digits and the results are rounded 
to the nearest 0.1 pm. Both the C++ and Maxima implementa- 
tions of the exa ct mapping and the test data are provided with 
GeographicLib (lKarnevil2010l) . 

The C++ implementation was checked by computing the 
maximum of the error in the forward mapping expressed as 
a true distance (i.e., dividing the error in the mapped space 
by the scale of the mapping) and the error in the reverse 
mapping (again expressed as a true distance). When imple- 
mented using double (resp. extended) precision, the maximum 
round-off error is 5r = 9nm (resp. 5 pm) over the whole 
ellipsoid (using the WGS84 parameters, a = 6 378137 m 
and / = 1/298.257 223 563). These ai-e consistent with 
6r ~ M/2P^^ indicating that the error is only about 8 times 
the limiting round-off error given in Sect. [T] The truncation 
error 5t, defined in Sect.|4] is zero for this method. 

Using the double precision implementation, the errors in 
the meridian convergence and scale at a particular point are 
bounded by 



. 1 

Skr 



2P-3 

1 



M ,/M\ 180° 



1 + 1.5 




where Sp and Sb are the geodesic distances from the point to 
the closest pole and closest branch point, respectively. These 
bounds were found empirically; however the form of the ex- 
pressions is determined by the nature of the singularities in 
the mapping. The term involving Sp arises because small er- 
rors in the position close to the pole may cause large changes 
in the convergence. Similarly the terms involving Sf, appear 
because of the singularity in the second derivative of the map- 
ping which causes the convergence and scale to vary rapidly 
near the branch point. 

The differences between my implementation and that of 
Do^ie£(1980) are as follows, (a) Dozier's starting guesses for 
Newton's method are based only on the limit e 0. New- 
ton's method then fails to converge in the neighborhood of the 
branch point (where ellipsoidal effects become large). In con- 
trast, 1 use different methods for computing the starting points 
in different regions which enables Newton's method to con- 
verge everywhere, (b) I modified several of the equations to 
improve the numerical accuracy. Without this, Dozier loses 
about half the precision in some regions, (c) 1 use published 
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FIG. 1 Graticules for the (a) spherical transverse Mercator, 
(b) Gauss-Kruger, and (c) Thompson projections. Here x and y are 
the easting and northing for the mappings. The eccentricity is e = 
(/ 1/199.5) and the mappings have been scaled so that the dis- 
tance from the equator to the north pole is unity. Thus A = 0° maps 
to the line a:: = and A = 90° maps to the line y = 1. The graticule 
is shown at multiples of 10° with 1° lines added in 80° < A < 90° 
andO < < 10°. 

algorithms to evaluate the special functions dBulirschl 1 19651; 
ICarlsonlll995h . (d) I compute the meridian convergence and 
scale, (e) Lastly, I provide an arbitrary precision implementa- 
tion (in Maxima) to allow the errors in the C-n- implementa- 
tion to be measured accurately. 

Figures [Ttb) and (c) show the graticule for the Gauss- 
Kriiger and Thompson projections. One eighth of the ellip- 
soid is shown in these figures, 0° < < 90° and 0° < A < 
90°, and (unlike the spherical transverse Mercator projection. 
Fig. [Tfa)) this maps to a finite area with the Gauss-Kriiger 
and Thompson projections. To obtain the graticule for the 
entire ellipsoid, reflect these figures in a; = 0, y = 0, and 
y = \. The eccentricity for these figures \s,e = and, in this 
case, the equator runs along y = until the branch point at 
A = Ao = 81° and then heads for y — 1; the point = 0°, 
A = 90° maps to finite points on y = 1. The Thompson 
mapping is not conformal at the branch point (where there's a 
kink in the equator), because dx/dw vanishes there. Similar 



figures are given bvlLei d 19761 Figs. 43^6). 

4. EXTENDING KRUGER'S SERIES 

There are several ways that the series for A, aj and /3j can 
be gene rated; here, I ad opt an approach which is close to that 
used bv iRriige^ (Il912l §5). (Alternative methods are to ex- 
pand Eqs. ( l30b and dSTI ) or to use the polar stereographic pro- 
jection instead of the Mercator projection as the starting point 
( IWallislll992l) .) In the limit 7] ^ rj' = (i.e., on the central 
meridian, A = 0), the quantities C, and C,' become the recti- 
fying and conformal latitudes respectively. (The conformal 
latitude (f)' was introduced in Sect. IH the rectifying latitude 
is linearly proportional to the distance along a meridian mea- 
sured from the equator.) The transformation between ( and C 
is thus given by the relation between the rectifying and con- 
formal latitudes extended to the complex plane. Thus ( is the 
meridian distance scaled to 7r/2, 

^(^) = 2^r (l^eW0)3/^ ^^' '''' 

where $ is the normal geographic latitude extended to the 
complex plane. The integral here can be expressed in terms 
of elliptic integrals as E{e) — E{Q, e), where E{(p, k) is 
the incomplete elliptic inte gral of the second kind with ar- 
gument (f) and modulus k (^DLMF, Eq. (19.2.5)), and Q = 
cot^^ ((1 — /) tan $) is the parametric co-latitude. Similarly, 
is merely Eq. (|4|i extended to the complex plane. 

C'($) = gd(gd~^ $ - etanh"^(esin$)). (33) 

The quantities and O are related to the Thompson projection 
variable w by 

$ ~ amw, O = &m{K(e) — w), 

(34) 

w = F($,e) = is:(e) -F(e,e), 

where am w is Jacobi's amplitude function ('dlmf', 5:22.16(1)1 
with modulus e and F{(t>, k) is the incomplete ellipt ic inte- 
gral of the first kind with argument cj) and modulus k (IdlmfI 
Eq. (19.2.4)). Substit uting E g. ( |34] i into Eqs. ^ and (O 
and using Eq. Q and lDLMH Eq. (22.16.31) gives Eqs. ( [30b 
and OU; this establishes t he equivale nce of Eqs. ( l32b and 
( l33b wi th the for mulation of'Le^ (Il976h . These equations are 
used bv lStuifberg en (2009) as the basis for an exact numerical 
method for the transverse Mercator proje ction. Th is is similar 
to (but rather simpler than) the method of lDozieJ (fl980) . 

The functions C,{^) and C'(^) analytic and so define 
conformal transformations. They can be expanded as a Tay- 
lor series in e^, or eq uivalently in n; I use the me thod of 
iLaerangd (Il770l §16) dWhittaker and WatsonL 1 19271 §7.32) 
to invert these series to give the inverse functions <I>(C) and 
<!>((■')■ For example, if 

C(<i>) = $ + 5('i>), 

where g{^) — 0{n), then the inverse function is 

$(C) = C + /i(C), 
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where 



MC) = E 



<I>=C 



Now compose C(*^'(C')) ^i^^ C'(^(C)) to provide the required 
series, Eqs. (|5]l and (|6]l. These manipulations were carried out 
using the algebraic tools provided by Maxima (2009). Little 
effort was expended to optimize this calculation since it only 
needs to be carried out once! (The expansion to order takes 
about 15 seconds.) At 8th order, the series for A is given 
by Eq. ( |29l ) and the series for aj, Eq. ( fT2] l. and /3j, Eq. ( fTTl ), 
become 

ai = in - ^n^ + A„3 , Ji 4 _ 127 5 , 7891 6 

^ ]_ r\ 1 1, rj/i/ V ^ G ' ^ lion"' '-)oo'*' I ■STonn 



16' 



180' 



288' 



37800 ' 



72161 7 _ 18975107 8 , 

T'^ r;r\or\onnri ~l 



387072 



50803200 



= l|n2 _ ^ _^ 2|1^5 _ 12||4|2n6 



1440 

13769 7 I 148003883 8 



630 



1935360 



28800 ' 



174182400 



_ ^„3 _ 103ti4 I 15061 5 , 167603 6 _ 67102379 7 



240 



140 



26880' 



181440 



29030400 ' 



79682431 8 
79833600 



161280 



168 



7257600 



49896 



40176129013 8 , 
7664025600 "• 



_ 34729 ^5 _ 3418889 ^6 _|_ 14644087 ^7 



80640 



1995840 



9123840 



2605413599 ^8 



' 622702080 
212378941^6 30705481^7 , 175214326799, 



319334400 



10378368 



58118860800 



_ 1522256789 7 _ 16759934899 8 



1383782400' 
_ 1424729850961 „ 



3113510400 



743921418240 



(35) 



and 



/Jl — 2"- 3"- 96" 360 

_ 5406467 7 , 7944359 8 , 

38707200 67737600 * ' ' ' ' 

^^2 T 45" 1440"" ' 105"' 

_i_ 51841 ^7 _|_ 24749483 ^8 _|_ . . . 



81 5 I 96199 6 
512 604800 



1118711 6 
3870720"' 



1209600 



348364800 



^3 = 



1L„3 _ 37_4 _ .209.„5 , 5569 „6 , 9261899 „7 

^ o/in'*' 4480 Qn'70n'*' ~r r^oncnonn'^ 



480 



840' 
6457463 8 
17740800 



90720 ' 



58060800 ' 



4397 4 1]_ 5 

161280 504 7257600 

I 324154477 „8 , 



830251 ^6 _j_ 466511 ^7 



2494800 



7664025600 ' 



4583_ 5 _ 108847 6 _ 8005831 7 



161280 



3991680' 



63866880' 



/36 = 



/37 = 
/38 



I 22894433 8 

124540416 ' 
20648693, ^6 



638668800 



16363163_^7 



518918400 



_ 2204645983 8 

12915302400 

219941297 7 
7 /t 



5535129600' 
191773887257 



497323811 8 
12454041600 



3719607091200'" '" ' ' ' ' 

Geographic Lib jKarnevlboiOl) includes the Maxima code for 
carrying out these expansions and the results of expanding the 
series much further, to order v?^ . 

Equations ( |29l ), ( [35] ), and ( [36l ) allow the Kxiiger method to 
be implemented to any order up to r? . In order to determine 



which order to use in a given application, it is useful to dis- 
tinguish the truncation error (the difference between the series 
evaluated exactly and the exact mapping) from the round-off 
error (the difference between the series evaluated at finite pre- 
cision and the series evaluated exactly). 

The truncation error was determined using Maxima's big- 
floats with a precision of 80 decimal digits. With Kriiger's 
series the error is principally a function of distance from the 
meridian (which is mainly a function of x) and depends only 
weakly on y. Defining s„j as the geodesic distance from 
the central meridian, I measure &t the maximum of the for- 
ward and reverse truncation errors (both expressed as true dis- 
tances) over all points with a given In Fig.|2l I plot 8t as a 
function of Sm, for truncations at various orders J (the small- 
est terms retained are n'). The errors rise monotonically with 
the distance from the central meridian. The branch point at 
(/) = 0° and A = Ao w 82.636° is at about = 9200 km. At 
this point the truncation error stops decreasing with increas- 
ing order indicating a lack of convergence in the series. (See 
Sect. |5]for a proof.) From a practical standpoint, the conver- 
gence is too slow to be useful for Sm ^ 8000 km. The trunca- 
tion errors in the meridian convergence and scale 5kt are 
well approximately by 



(57t = 2 J sec(sm/a) 



(St 180° 



a TT 
<St 



— — — 2J sec{s„i/a) — . 
k a 

Here, the factor 2 J arises from the differentiation performed 
to give Eqs. ( |23] l and ( |26] l and the term sec(s„i/a) is the scale 
(in the spherical limit) necessary to convert the errors in posi- 
tion to errors in ( or (^'. 

Figure |2] epitomizes the advantages of Kriiger's over 
Thomas' series. The equivalent figure for truncation error for 
the latter series would use the longitude relative the central 
meridian for the abscissa instead of the distance from the cen- 
tral meridian. At high latitudes, the longitude difference be- 
comes large even for modest distances from the central merid- 
ian; this explains the large errors in the results from Geotrans 
in the Greenland example at the end of Sect.|2l 

Round-off errors need to be considered also when imple- 
menting the method with floating-point arithmetic. Evaluating 
the mapping using the formulas given in Sect.|2]adds 4.2 nm 
or L9pm to the truncation error depending on the precision 
(see the dashed lines in Fig. These round-off errors may 
be expressed as Sr ~ A//2p^^; i.e., they are about half of 
those for the exact algorithm, because the series method in- 
volves fewer operations. Thus with double (resp. extended) 
precision and the series truncated at J = 6 (resp. J — 8) 
the overall error is less than 5 nm (resp. 2 pm) provided that 
Sm < 3900 km (resp. 4200 km). The C-n- implementation 
of the mapping based on the extended Kriiger series (taken to 
order n^) is included in GeographicLib (Karnev, 2010). 

At greater distances from the central meridian, truncation 
errors become large; thus the 6th order series has an error of 
about 1 mm at s„i — 7600 km (see Fig. |2|i. The truncation 
error can be decreased by increasing the number of terms re- 
tained in the series. However, if high accuracy is required for 
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FIG. 2 Truncation and round-off errors, 5t and Sr, in the Kriiger 
series for the transverse Mercator projection as a function of the dis- 
tance from the central meridian Sm- The series is truncated at various 
orders J from 2 to 12. The solid lines show the St. Also shown in 
dashed lines are the combined round-off and truncation errors, St+S,. 
when the algorithm is implemented with floating-point numbers with 
p bits of fraction. For p = 53 (double) and p = 64 (extended), the 
truncation levels are set to J = 6 and J — S respectively. The 
WGS84 ellipsoid is used. 

Sm ^ 4000 km, it's probably safer to use the exact algorithm 
whose error is less than 9 nm for the whole spheroid. 

The round-off errors in the meridian convergence and scale 
are bounded by 

these should be added to Sjt and Skf. These expressions have 
the same form as those for the exact algorithm, except that the 
terms involving Sb are omitted because the truncation eiTor 
dominates nea r the branch point. 

Previously, Engsager and Pode3 (l2007h extended Kriiger's 
series to 7th order. However they give a less rigorous esti- 
mate of the error. In addition, they give separate series for 
the meridian convergence and scale, while I advocate follow- 
ing Kriiger's simpler prescription of differentiating the same 
series used to carry out the mapping. They also give a series 
expansion for the transformation between latitude and con- 
formal latitude which leads to a somewhat faster code (see 
Sect. |6|i; however, I prefer the simplicity of evaluating and 
inverting Eq. (|7]i directly. 



5. PROPERTIES FAR FROM THE CENTRAL MERIDIAN 

In this section, I explore the behavior of the mapping in the 
vicinity of the branch point at = 0° and A = Aq. First it 
should be remarked that this is far from the central meridian 
and that the scale there is fco / e; so this is not in the domain 
where the mapping is very useful. Nevertheless, it is benefi- 
cial to have a complete understanding of a mapping; for exam- 
ple, this was necessary in making the exact algorithm robust. 
Here I d escribe the mapping in terms of the Tho mpson map- 
ping, w (lLeg.ll976l §§54-55). ' Konig and Weis3 (l951) offer 
a complementary picture based on the complex latitude $. 

The transverse Mercator projection is defined by its prop- 
erties on the central meridian and the condition of confor- 
mality. The mapping i s defin ed by analytically continuing 
jWhittaker and Watsonl 1 19271 Chap. 5) the mapping away 
from the central meridian. The process continues until the 
branch point is encountered (where the condition of analyt- 
icity fails). The branch point coiTesponds to x = Xo = * Aq, 
C = Co = iiK'-E')Tr/{2E),andw ^ wq = iK',whsTsE = 
E{e), E' = E{Vl-e^), K = K{e), K' = K{Vl~S^), 
and K{k) is the com plete elliptic integral of the first kind with 
modulus k (IdlmfI Eq. (19.2.8)). The lowest order terms in 
the expansions of x and C, about the branch point are 

X = - e") [w^ + ^(1 + e2)^5 + ...], 

C - - e^) [w^ + 1(2 - e^)w' + • • • ] V(2i?), 

where x — X~ Xo^ C = C, ~ Co^ and w — w — wq- Efiminating 
w from these equations gives 

The value of C, will depend on how the | power is taken. Pick- 
ing the complex phase of ix in the interval (— tt, tt] gives the 
principal value. Equivalently, the value of C(x) can be made 
single valued by placing "cuts" on the equator in the longi- 
tude ranges (1 - e)90° < |A| < (1 + e)90° which act as 
impassable barriers during the process of analytic continua- 
tion. This represents the "standard" convention for mapping a 
geographic position to the Gauss-Kriiger projection since the 
sign of the northing matches the sign of the latitude (with the 
equator mapping to non-negative northings). This convention 
coiTesponds to Fig.IIlb) (after suitab le reflections to cover the 
ellipsoid), to Lee ( 1976, Fi g. 46) , to iKonig and Weisd (Il95l[ 
Fig. 55(b)), and to .Ludwigi(ll943i p. 214). 

From the form of the mapping near the branch point, it is 
clear that the Kriiger series does not converge for (p — and 
A > Aq because, from Eq. ( |37] ). ( is complex under these 
conditions but all the terms in the series, Eq. (|5]), are pure 
imaginary. Delineating the precise boundary for convergence 
of the series and its inverse, Eq. (|6]l, requires an analysis of the 
problem for complex e and is beyond the scope of this work. 

It is possible to extend the mapping by moving to the 
"right" of the equator in Fig. [Tfb)- If the complex phase of 
ix includes the interval [0, |7r], an "extended" domain for the 
mapping may be defined by the union of 0° < (j> < 90°, 
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FIG. 3 Extended transverse Mercator projection, (a) shows the 
graticule and (b) the convergence and scale for Gauss-Kriiger pro- 
jection, (c) shows the graticule for the Thompson projection. The 
ellipsoid parameters are the same as for Fig.[T] The graticules in (a) 
and (c) are the same as in Fig.[TJb) with the addition of 1° lines for 
80° < A < 90° and -10° < </) < 0°. In (b), the lines emanating 
from the top left comer are lines of constant meridian convergence, 7, 
at 10° intervals. The dog-legged line joining (0, 1), (0, 0), (1.71, 0), 
and (1.71, —00) represents 7 = 0°. The line y — 1 gives 7 = 90°. 
The other lines (running primarily vertically in the figure) are lines 
of constant scale k. The solid lines show integer values of k for 
1 < fc < 15 and multiples of 5 for 15 < < 35. The line segment 
joining (0, 0) and (0, 1) gives fc = 1. The dashed lines show lines of 
constant k at intervals of 0.1 for 1 < fc < 2. 

0° < A < 90° and -90° < </) < 0°, Ao < A < 90°. 
The rule for analytic continuation is that the second region is 
reached by a path from the central meridian which goes north 
of the branch point. This is equivalent to placing the cut so 
that it emanates from the branch point in a south-westerly di- 
rection. Following this prescription, the range of the mapping 
now consists of the union of < x, < y/{kQa) < E and 
K' - E' <x/{kaa),y <Q. 

Figures [3ja) and (b) illustrate the properties of the Gauss- 
Kriiger projection in this extended domain. These figures use 
an ellipsoid with eccentricity e — (as in Fig. [1) and with 
a = 1/E = 0.6382 and fco = 1. The branch point then lies 



at = 0°, A = 81° orx= [K' - E')/E w 1.71, y = 0. 
Symmetries can now be employed to extend the mapping with 
arbitrary rules for how to circumvent the branch point. The 
symmetries are equivalent to placing mirrors on the four lines 
segments: 0<x, ?/ = l;a; = 0, 0<?/<l;0<a;< 
1.71, ?/ = 0; and x = 1 .71, y < 0. Compare Fig.[3la) with 
'Konig and WeisS ([1951', Fig. 53(b)). 

Figure[3c) shows the graticule of the Thompson projection 
in the extended domain; the range of this mapping is the rect- 
angular region shown, < x/{kQa) < K', < y/{kQa) < 
K. The extended Thompson projection has reflection sym- 
metry on all the four sides of Fig. Oc). In transforming from 
Thompson to Gauss-Kriiger, the right angle at the lower right 
corner of Fig. Oc) expands by a factor of 3 to 270° to pro- 
duce the outside corner at x — 1.71, y = Q in Fig.[3ta). The 
top right corner of Fig. Oc) represents the south pole and this 
is transformed to infinity in the extended Gauss-Kriiger pro- 
jection. Despite the apparent similarities, the behavior of the 
extended Thompson projection near the north and south poles 
(the top left and top right corners in Fig. Oc)) is rather dif- 
ferent. Although the mapping at north pole is conformal, the 
mapping at the south pole in the extended domain is not. The 
difference in longitude between the two meridians represent 
by the top and right edges of Fig. [3c) is 90°e instead of 90°. 

My implementation of the exact mapping provides the op- 
tion of using the extended domain. The round-off eiTors 
quoted in Sect.[3](9 mii for double precision and 5 pm for ex- 
tended precision) apply to the extended domain for </) > —15°. 
Beyond this line, the errors grow because of the contraction of 
w space near the south pole (at (j> = —58°, the error is about 
1 mm). 

6. CONCLUSION 

The algorithms presented here allow the transverse Mer- 
cator projection to be computed with an accuracy of a few 
nanometers. Implement ations of these algorithms are in- 
cluded in GeographicLib dKarnevlEonl which also provides 

(a) the set of test data used to check the implementations, 

(b) Maxima code for the exact mapping (with arbitrary preci- 
sion), (c) Maxima code for generating the Kriiger series to ar- 
bitrary order, and (d) the Kriiger series to 30th order. The web 
page http://geographiclib.sf.net/tm.html provides quick links 
to all these resources. 

The work described in this paper mad e heavy use of the 
computer algebra system iMaxim i (120091) both for carrying 
out the series expansions for Kriiger's method and for gen- 
erating the high accuracy test data. The latter is invaluable 
when developing complex algorithms with an accuracy close 
to machine precision. Other computer algebra systems offer 
similar capabilities; but Maxima is one of the few that is free. 

My emphasis in developing these algorithms was in their 
accuracy. Nevertheless the resulting implementations are rea- 
sonably fast. On a 2.66 GHz Intel processor and compiled 
with g++, the time for the mappings implemented with the 6th 
order series method is 1.91 /is; this is the combined time for a 
forward and a reverse mapping including the computation of 
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the convergence and scale in each case. This time is insen- 
sitive to number of terms retained in the sum due to the effi- 
ciency of Clenshaw summation — changing this to 4 (resp. 8) 
decreases (resp. increases) the time by only 1%. Skipping the 
calculation of the convergence and scale reduces the time by 
15%. Using a trigonometric series and Clenshaw summation 
for the conversions between geographic and conformal lati- 
tude (as proposed by Engsager and Poder (2007)) decreases 
the time by 18%. The exact algorithms (which are accurate 
over the entire ellipsoid) are 5-6 times slower. The 6th order 
series method is comparable in speed to Geotrans 3.0 even 
though the latter is much less accurate and does not return the 
convergence and scale. 

Here are some recommendations for users of the transverse 
Mercator projec tion. Do not u se algorithms based on the for- 
mulas given by iThomasI ( Il952) — they are unnecessarily in- 
accurate. Instead use the Kriiger series, truncating Eqs. ( [29] ). 
( [35] l. and ( [36l ) to order n^. With double precision, this gives an 
accuracy of 5 nm for distances up to 3900 km from the cen- 
tral meridian. If the mapping is needed at greater distances 
from the central meridian, use the algorithm based on the ex- 
act mapping (an accuracy of 9nm). If greater accuracy is 
needed, use extended precision with either method (extend- 
ing the series method to 8th order). When implementing these 
algorithms, use the test set to verify that the errors are compa- 
rable with those given here. 
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